#include <vector>
#include <iostream>
#include <iomanip>
#include <cmath>
#include "RI.h"
#include "Grid.h"
#include "Solve.h"

using namespace std;

class pde_u : public Function {
public:
    vector<double> operator()(vector<double> x) {
        vector<double> result(x.size(),0.0);
	for (int i = 0; i < x.size(); i++)
	{
            result[i] = -sin(x[i])*exp(sin(x[i])) + cos(x[i])*cos(x[i])*exp(sin(x[i]));
        }
        return result;
    }
};

class du : public Function {
public:
    vector<double> operator()(vector<double> x) {
        vector<double> result(x.size(),0.0);
	for (int i = 0; i < x.size(); i++)
	{
            result[i] = cos(x[i])*exp(sin(x[i]));
        }
        return result;
    }
};

class u : public Function {
public:
    vector<double> operator()(vector<double> x) {
        vector<double> result(x.size(),0.0);
	for (int i = 0; i < x.size(); i++)
	{
            result[i] = exp(sin(x[i]));
        }
        return result;
    }
};

int main() {
    //Grid grid1d(1, 5);      // 创建一维网格，分隔数为5
    //grid1d.display();
    
    //Grid grid2d(2, 3);      // 创建二维网格，分隔数为3
    //grid2d.display();
    //Grid grid3d(3, 4);      // 创建三维网格，分隔数为4，非法维度
    MatrixVectorOperator MVO;
    int N = 32;
    pde_u f;
    u u;
    cout << endl;
    int k = 0;
    vector<vector<double>> r1(4), r2(4);
    for (N = 32; N <= 256; N = N*2)
    {
	Grid grid1D(1, N); //1维N区间网格
	vector<double> x = grid1D.grid1D;
	vector<double> U0 = u(x), U(N-1), F0 = f(x), F(N-1),v0(N-1),v1(N-1),v2(N-1);
	for (int i = 1; i < N; i++)
	{
	    U[i - 1] = U0[i];
	}
	Solver S2(N,f,v0,"D","FullWeight","Linear","FMG");
	Solver S1(N,f,v0,"D","FullWeight","Linear","V-cycle");
	v1 = S1.solve(U0[0],U0[N],20, 1e-8);
	v2 = S2.solve(U0[0],U0[N],20, 1e-8);
	r1[k] = S1.Re;
	r2[k] = S2.Re;
	k ++;
    }
    cout << "V-cycle 迭代时残差的无穷范数取log_2" << endl;
    cout << "iter" << "\t" << "N = 32" << "\t" << "N = 64" << "\t" << "N = 128" << "\t" << "N = 256" << endl;
    for (int i = 0; i < 9; i++)
    {
    	cout << i ;
	for (int j = 0; j < 4; j++)
    	{
	    //cout << "\t" << log(r1[j][i])/log(2);
	    cout << "\t" << fixed << setprecision(3) << log(r1[j][i])/log(2);
    	}
	cout << endl;
    }
    cout << endl;
    cout << "FMG 迭代时残差的无穷范数取log_2" << endl;
    cout << "iter" << "\t" << "N = 32" << "\t" << "N = 64" << "\t" << "N = 128" << "\t" << "N = 256" << endl;
    for (int i = 0; i < 9; i++)
    {
    	cout << i ;
	for (int j = 0; j < 4; j++)
    	{
	    //cout << "\t" << log(r1[j][i])/log(2);
	    cout << "\t" << fixed << setprecision(3) << log(r2[j][i])/log(2);
    	}
	cout << endl;
    }
    return 0;
}
